// 版权所有2010 Go作者。保留所有权利。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// 下面的原始C代码、长注释和常量
// 来自http:
// go代码是原始C.
// 
// tgamma的简化版本。c 
// 
// Gamma函数
// 
// ；
// extern int SINGGAM；
// 
// y=tgamma（x）；chuang jian defg
// 符号正确，符号（+1或-1）也是
// 在名为signgam的全局（外部）变量中返回。
// 此变量也由对数伽马函数填写
// 函数lgamma（）。
// 
// 参数| x |<=34通过递归减少，函数
// 在
// 区间（2,3）中近似为6/7阶有理函数。大型争论由斯特林的
// /公式处理。使用
// 一个反射公式将大的否定论点变为正。
// 
// 准确度：
// 
// 相对误差：
// 算术域#试验峰值均方根值
// DEC-34，34 10000 1.3e-16 2.5e-17 
// IEEE-170，-33 20000 2.3e-15 3.3e-16 
// IEEE-33，33 20000 9.4e-16 2.2e-16 
// IEEE 33、171.6 20000 2.3e-15 3.2e-16 
// 
// 由于指数函数的误差放大，测试范围外参数的误差将更大
// 。
// 
// Cephes数学库发行版2.8:2000年6月
// 版权所有1984、1987、1989、1992、2000作者Stephen L.Moshier 
// 
// 自述文件http:
// 本档案中的一些软件可能来自《数学函数的方法和程序》（Prentice Hall或Simon&Schuster 
// International，1989）一书，也可能来自Cephes数学图书馆，一本
// 商业产品。无论哪种情况，它都是作者的版权。
// 您在这里看到的内容可以免费使用，但不提供任何支持或
// 保证。
// 
// 本书中两个已知的印刷错误在这里的
// 伽马函数和不完整的beta 
// 积分的源代码列表中修复。
// 
// Stephen L.Moshier 
// moshier@na-网。奥尔。gov 

var _gamP = [...]float64{
	1.60119522476751861407e-04,
	1.19135147006586384913e-03,
	1.04213797561761569935e-02,
	4.76367800457137231464e-02,
	2.07448227648435975150e-01,
	4.94214826801497100753e-01,
	9.99999999999999996796e-01,
}
var _gamQ = [...]float64{
	-2.31581873324120129819e-05,
	5.39605580493303397842e-04,
	-4.45641913851797240494e-03,
	1.18139785222060435552e-02,
	3.58236398605498653373e-02,
	-2.34591795718243348568e-01,
	7.14304917030273074085e-02,
	1.00000000000000000320e+00,
}
var _gamS = [...]float64{
	7.87311395793093628397e-04,
	-2.29549961613378126380e-04,
	-2.68132617805781232825e-03,
	3.47222221605458667310e-03,
	8.33333333333482257126e-02,
}

// 用斯特林公式计算的伽马函数。
// 这对结果必须相乘才能得到实际答案。
// 乘法留给调用者，因此，如果小心，调用者可以避免172<=x<=180的无穷大。
// 多项式对33<=x<=172有效；较大的值仅用于倒数运算
// 并产生非规范化浮点。精度较低的
// 掩盖了多项式中的任何不精确性。
func stirling(x float64) (float64, float64) {
	if x > 200 {
		return Inf(1), 1
	}
	const (
		SqrtTwoPi   = 2.506628274631000502417
		MaxStirling = 143.01608
	)
	w := 1 / x
	w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
	y1 := Exp(x)
	y2 := 1.0
	if x > MaxStirling { // avoid Pow（）overflow 
		v := Pow(x, 0.5*x-0.25)
		y1, y2 = v, v/y1
	} else {
		y1 = Pow(x, x-0.5) / y1
	}
	return y1, SqrtTwoPi * w * y2
}

// Gamma返回x的Gamma函数。
// 
// NaN 
func Gamma(x float64) float64 {
	const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620 
	// 特殊情况
	switch {
	case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
		return NaN()
	case IsInf(x, 1):
		return Inf(1)
	case x == 0:
		if Signbit(x) {
			return Inf(-1)
		}
		return Inf(1)
	}
	q := Abs(x)
	p := Floor(q)
	if q > 33 {
		if x >= 0 {
			y1, y2 := stirling(x)
			return y1 * y2
		}
		// 注意：x是负数，但（上面选中）不是负整数，
		// 因此x必须足够小，才能在转换为int64的范围内。
		// 如果| x | were>=2⁶它必须是一个整数。
		signgam := 1
		if ip := int64(p); ip&1 == 0 {
			signgam = -1
		}
		z := q - p
		if z > 0.5 {
			p = p + 1
			z = q - p
		}
		z = q * Sin(Pi*z)
		if z == 0 {
			return Inf(signgam)
		}
		sq1, sq2 := stirling(q)
		absz := Abs(z)
		d := absz * sq1 * sq2
		if IsInf(d, 0) {
			z = Pi / absz / sq1 / sq2
		} else {
			z = Pi / d
		}
		return float64(signgam) * z
	}

	// 减少参数
	z := 1.0
	for x >= 3 {
		x = x - 1
		z = z * x
	}
	for x < 0 {
		if x > -1e-09 {
			goto small
		}
		z = z / x
		x = x + 1
	}
	for x < 2 {
		if x < 1e-09 {
			goto small
		}
		z = z / x
		x = x + 1
	}

	if x == 2 {
		return z
	}

	x = x - 2
	p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
	q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
	return z * p / q

small:
	if x == 0 {
		return Inf(1)
	}
	return z / ((1 + Euler*x) * x)
}

func isNegInt(x float64) bool {
	if x < 0 {
		_, xf := Modf(x)
		return xf == 0
	}
	return false
}
